
/*
 * TravelingSantaPaths.c
 *
 * Traveling Santa Problem - Paths chooser within a region.
 *
 *  Written by "ACME Labs" team (Anton Sletner & Todd Fisher)
 *   ("ACME Labs" as used here is only a kaggle team name.)
 *   for  kaggle  Traveling Santa Problem  contest on
 *   kaggle  ( www.kaggle.com )  ( 20121214 - 20130118 )
 *   ( Predictive Data Modeling Crowdsourcing )
 *
 *  If contest winner,
 *   Released under the GNU GPL v3 Open Source License:
 *   http://www.gnu.org/licenses/gpl.html
 *
 *  Like the Traveling Salesman Problem except that Santa needs
 *   two disjoint paths to all cities/nodes instead of one, and
 *   with no requirement to return to any starting city/node.
 *
 *  Usually contests at kaggle are predictive data modeling
 *   problems, but this one is an optimization problem instead.
 *
 *  Started  20121229  ALS
 *
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include "santas-tour.h"

void print_basic_nodes( basic_node *list )
{
    basic_node *nd;
    printf( "      #  |  number  |         x         y\n" );
    printf( "-----------------------------------------\n" );
    int c = 1;
    for( nd = list; ; nd++, c++ )
    {
        printf( "%7d  | ", c );
        printf( "%7d  |  ", nd->number );
        printf( "%8.2lf, %8.2lf \n", nd->x, nd->y );

        if( NULL == nd->next )  break;
    }
}

static double minx, maxx, miny, maxy, xdiff, ydiff;
static int nodecnt, xcnt, ycnt, maxcnt;

static int GREY;
static extra *extras;
static grey_node *greygrid;
static int minvirtx, maxvirtx, minvirty, maxvirty;
static int mingreyx, maxgreyx, mingreyy, maxgreyy, greyxcnt, greyycnt;

static int QTP2, QTDEPTH;
static int xqtadj, yqtadj;  // Centering stuff -- adj can be neg.
static int xqtoff, yqtoff;

void setMinMaxXY( basic_node *on_x, basic_node *on_y )
{
    minx = on_x->x;
    miny = on_y->y;

    int c=0;
    for( ; on_x[c].next != NULL; c++ );
    maxx = on_x[c].x;

    nodecnt = c+1;

    maxy = on_y[c].y;

    xdiff = maxx -minx;
    ydiff = maxy -miny;

    xcnt  = xdiff +1.0;
    ycnt  = ydiff +1.0;

    maxcnt = (xcnt>=ycnt)*xcnt + (ycnt>xcnt)*ycnt; 

    if( NULL == (extras = (extra *)malloc( sizeof(extra)*nodecnt )) )
    {
        fprintf( stderr, "Memory overflow on extras creation.\n" );
        exit(0);
    }
    for( c=0; c<nodecnt; c++ )  extras[c].number = c;

    // Extend virtual bounds to match coarser GREY bounds if necessary.
    GREY = (maxcnt * 0.005 +0.5);
    minvirtx = (((int)minx) / GREY) * GREY;
    maxvirtx = ((int)maxx) +(GREY -(((int)maxx) % GREY)) -1;
    minvirty = (((int)miny) / GREY) * GREY;
    maxvirty = ((int)maxy) +(GREY -(((int)maxy) % GREY)) -1;

    // Set coarse GREY block vals.
    mingreyx = minvirtx / GREY;
    maxgreyx = maxvirtx / GREY;
    mingreyy = minvirty / GREY;
    maxgreyy = maxvirty / GREY;

    greyxcnt = maxgreyx -mingreyx +1;
    greyycnt = maxgreyy -mingreyy +1;

    // To create virtual coarse grey grid rectangle overlay of "grey nodes".
    if( NULL == (greygrid = (grey_node *)malloc( sizeof(grey_node)*greyxcnt*greyycnt )) )
    {
        fprintf( stderr, "Memory overflow on grey grid alloc.\n" );
        exit(0);
    }


    // (Extend) x and y to the next power of 2 of the larger.
    //  (All nodes are assumed to be at integer coords, so we (extend) to a pow of 2
    //   side length square, so our quadtree doesn't have to deal with any overlaps.)
    QTP2 = 1;
    int sh = maxcnt, nz = 1;
    for( c=0; sh!=0 ; c++ ){  if( (sh&1) && (sh>>1) ) nz=0;  sh>>=1;  QTP2<<=1;  }
    QTP2 >>= nz;  // Shift down 1 if already a pow of 2.

    QTDEPTH = 0;
    int t;
    for( t=QTP2; t>1; t>>=1 ) QTDEPTH++;

    xqtoff = (QTP2>>1) -(xcnt>>1);
    yqtoff = (QTP2>>1) -(ycnt>>1);
    xqtadj = xqtoff -minx;
    yqtadj = yqtoff -miny;

    fprintf( stderr, "Min X: %lf, Max X: %lf, Min Y: %lf, Max Y: %lf  Nodes: %d\n", minx, maxx, miny, maxy, nodecnt );
    fprintf( stderr, "MinvX: %d, MaxvX: %d, MinvY: %d, MaxvY: %d  GREYbox: %d \n", minvirtx, maxvirtx, minvirty, maxvirty, GREY );
    fprintf( stderr, "QadjX: %d, QadjY: %d, QoffX: %d, QoffY: %d  QTP2: %d  QTDEPTH: %d \n", xqtadj, yqtadj, xqtoff, yqtoff, QTP2, QTDEPTH );
}

// Rasterize as rectangle
basic_node *scanAlongX( basic_node *on_y )
{
    #define TAIL 100
    double avggap, x, y, yi, pix[TAIL], initfill;
    int    gapi, bix = 0, xpix = minvirtx /* minx */, first = 1;
    int    gline, gy, gxpix = 0, gxcnt = GREY;

    basic_node *px = on_y, *next = px->next;

    // Treat as a rectangle picture on init.
    x  = px->x;
    y  = px->y;
    yi = y;
    gy    = yi / GREY;
    gline = (gy -mingreyy) * greyxcnt;
    gapi      = x -minvirtx /* minx */;  // Pre-scan count from edge of picture to first node.
    initfill  = ((NULL != next) && (next->y == y)) ? 1.0/(next->x -x) : 1.0/xcnt;
    int c, d, end = 0;
    for( c=TAIL-1; c>=0; c-- )  pix[c] = 0.0;

    while( y == yi )
    {
        for( c=gapi; c>=end ; c-- )  // Treat any gap as pixel(s).
        {
            if( 0 == c )  pix[bix] = 1.0;  // Note next node.
              else        pix[bix] = 0.0;  // Note gap to next node.

            // Running average gap.
            for( d=TAIL-1, avggap=0.0; d>=0; d-- )  avggap +=pix[d];
            avggap /= ((gapi>TAIL) ? gapi+1 : TAIL);

            if( 0 == c )
            {
                if( first ){  avggap = initfill;  pix[bix] = 0.0;  first=0;  }

                extras[px->number].grey_x = avggap;
//                fputc( 'X', stderr );
//              if( avggap == 1.0 )
//                printf( "y: %lf  xpix: %d:  avg: %lf\n", y, xpix, avggap );

                greygrid[gline+gxpix].nodecnt++;
            }

            
//            printf( "y: %lf  xpix: %d:  avg: %lf\n", y, xpix, avggap );

            // Shift
            bix = (bix +1) % TAIL;

            xpix++;

            gxcnt--;
            if( 0 == gxcnt ){

                if( ((int)yi % GREY) == (GREY -1) )  // printf( "gx:%d gy:%d:  cnt: %d\n", gxpix +mingreyx, gy, greygrid[gline+gxpix].nodecnt );

                gxpix++;  gxcnt = GREY;  }
        }

        next = px->next;

        if( (NULL == next) || (next->y != y) )  // Scan past last node to edge of picture.
            if( (gapi = maxvirtx /* maxx */ -(xpix -1)) != 0 )
            {
                end = 1;
                continue;
            }

        if( NULL == next )  break;

        gapi = next->x -x -1.0;
        px++;
        x = px->x;
        y = px->y;
    }

    if( NULL == px->next )  return  NULL;
    return  px;
}

// Rasterize as rectangle
basic_node *scanAlongY( basic_node *on_x )
{
    #define TAIL 100
    double avggap, y, x, xi, pix[TAIL], initfill;
    int    gapi, bix = 0, ypix = minvirty /* miny */, first = 1;
    int    gline, gx, gypix = 0, gycnt = GREY;

    basic_node *py = on_x, *next = py->next;

    // Treat as a rectangle picture on init.
    y  = py->y;
    x  = py->x;
    xi = x;
    gx = xi / GREY;
    gline = (gx -mingreyx) * greyycnt;
    gapi     = y -minvirty /* miny */;  // Pre-scan count from edge of picture to first node.
    initfill = ((NULL != next) && (next->x == x)) ? 1.0/(next->y -y) : 1.0/ycnt; 
    int c, d, end = 0;
    for( c=TAIL-1; c>=0; c-- )  pix[c] = 0.0;

    while( x == xi )
    {
        for( c=gapi; c>=end ; c-- )  // Treat any gap as pixel(s).
        {
            if( 0 == c )  pix[bix] = 1.0;  // Note next node.
              else        pix[bix] = 0.0;  // Note gap to next node.

            // Running average gap.
            for( d=TAIL-1, avggap=0.0; d>=0; d-- )  avggap +=pix[d];
            avggap /= ((gapi>TAIL) ? gapi+1 : TAIL);

            if( 0 == c )
            {
                if( first ){  avggap = initfill;  pix[bix] = 0.0;  first=0;  }

                extras[py->number].grey_y = avggap;
//                fputc( 'Y', stderr );
//              if( avggap == 1.0 )
//                printf( "x: %lf  ypix: %d:  avg: %lf\n", x, ypix, avggap );
            }

//            printf( "x: %lf  ypix: %d:  avg: %lf\n", x, ypix, avggap );

            // Shift
            bix = (bix +1) % TAIL;

            ypix++;

            gycnt--;
            if( 0 == gycnt ){  gypix++;  gycnt = GREY;  }

        }

        next = py->next;

        if( (NULL == next) || (next->x != x) )  // Scan past last node to edge of picture.
            if( (gapi = maxvirty /* maxy */ -(ypix -1)) != 0 )
            {
                end = 1;
                continue;
            }

        if( NULL == next )  break;

        gapi = next->y -y -1.0;
        py++;
        y = py->y;
        x = py->x;
    }

    if( NULL == py->next )  return  NULL;
    return  py;
}

void find_paths( basic_node *sorted_on_x, basic_node *sorted_on_y, node **path1, node **path2 )
{
//    printf( "MOO!\n" );

//    print_basic_nodes( sorted_on_y );

    setMinMaxXY( sorted_on_x, sorted_on_y );

// --
/*
    fprintf( stderr, "Creating grey grid. \n" );

    basic_node *next_y = sorted_on_y, *next_x = sorted_on_x;

    do
    {
        next_y = scanAlongX( next_y );

    } while( NULL != next_y );

    do
    {
        next_x = scanAlongY( next_x );

    } while( NULL != next_x );
*/
// --

    fprintf( stderr, "Creating quadtree. \n" );

    qt_box box;
    qt_box_init( &box, 0.0, 0.0, QTP2, QTP2 );
    quadtree *qt = qt_create( box, nodecnt, QTDEPTH /*, overlap */ );

    basic_node *np = sorted_on_x;
    int c;
    for( c=1; ; np++, c++ )
    {
        // fprintf( stderr, "Node: %d \n", np->number );
        qt_insert( qt, np, &box );

        if( NULL == np->next )  break;
    }
    fprintf( stderr, "Inserted %d.  Done creating quadtree -- Now destroying. \n", c );

    qt_destroy( qt );

    fprintf( stderr, "Done destroying quadtree. \n" );
}

// --- Quadtree stuff -- Unfinished.
//      For now we'll only use on a 2**n by 2**n grid, eliminating the
//       need for dealing with overlaps--All splits are to integer widths.

int qt_box_is_valid( qt_box *qb )
{
    return  ((qb->xmin < qb->xmax) && (qb->ymin < qb->ymax)) ? 1 : 0;
}

int qt_box_is_inside( const qt_box *first, const qt_box *second )
{
    return  ((second->xmin < first->xmin) && (second->xmax > first->xmax) &&
             (second->ymin < first->ymin) && (second->ymax > first->ymax)) ? 1 : 0;
}

int qt_box_overlaps( const qt_box *first, const qt_box *second )
{
    return  ((first->xmin > second->xmax) || (first->xmax < second->xmin) ||
             (first->ymin > second->ymax) || (first->ymax < second->ymin)) ? 0 : 1;
}

double qt_box_width( const qt_box* qb )
{
    return  (qb->xmax - qb->xmin);
}

double qt_box_height( const qt_box* qb )
{
    return  (qb->ymax - qb->ymin);
}

void qt_box_init( qt_box *qb, double xmin, double ymin, double xmax, double ymax )
{
    qb->xmin = xmin; qb->ymin = ymin; qb->xmax = xmax; qb->ymax = ymax;
    assert( qt_box_is_valid( qb ) );
}

void qt_box_inflate( qt_box *qb, double dx, double dy )
{
    qb->xmin -= dx/2.0;
    qb->xmax += dx/2.0;
    qb->ymin -= dy/2.0;
    qb->ymax += dy/2.0;
}

qt_node *qt_node_create( const qt_box *box )
{
    qt_node *node = NULL;

    if( NULL == (node = (qt_node*)malloc( sizeof(qt_node) )) )
    {
        fprintf( stderr, "Memory overflow on quadtree node creation.\n" );
        exit(0);
    }

    memset( node->child, 0, sizeof(node->child) );
    memcpy( &(node->box), box, sizeof(qt_box) );
    node->pconflicts = NULL;

    return  node;
}

quadtree *qt_create( qt_box box, int nodecnt, int depth /*, float overlap */ )
{
    quadtree* qt = NULL;

    if( NULL == (qt = (quadtree*)malloc( sizeof(quadtree) )) )
    {
        fprintf( stderr, "Memory overflow on quadtree init of creation.\n" );
        exit(0);
    }

    qt->depth   = depth;
 /* qt->overlap = overlap;

    qt_box_inflate( &box, qt_box_width(&box)*overlap, qt_box_height(&box)*overlap );
 */
    qt->root = qt_node_create( &box );

    if( !qt->root )
    {
        free( qt );
        qt = NULL;
    }

    return  qt;
}

void qt_node_destroy( qt_node* node )
{
    if( node->child[NE] != NULL )
    {
        qt_node_destroy( node->child[NE] );
        qt_node_destroy( node->child[NW] );
        qt_node_destroy( node->child[SW] );
        qt_node_destroy( node->child[SE] );
    }

    free( node );
}

void qt_destroy( quadtree *qt )
{
    if( qt && qt->root )
    {
        qt_node_destroy( qt->root );
        free( qt );
    }
}

void qt_box_split( const qt_box *qb, qt_box *ne, qt_box *nw, qt_box *se, qt_box *sw /*, float overlap */ )
{
    double dx = qt_box_width(qb)  * (1.0 /*+ overlap */ )/2.0;
    double dy = qt_box_height(qb) * (1.0 /*+ overlap */ )/2.0;

 /* assert( overlap >= QBOX_OVERLAP_MIN-0.0001 ); */
 /* assert( overlap <= QBOX_OVERLAP_MAX+0.0001 ); */

    qt_box_init( ne, qb->xmax-dx, qb->ymax-dy, qb->xmax,    qb->ymax );
    qt_box_init( nw, qb->xmin,    qb->ymax-dy, qb->xmin+dx, qb->ymax );    
    qt_box_init( sw, qb->xmin,    qb->ymin,    qb->xmin+dx, qb->ymin+dy );
    qt_box_init( se, qb->xmax-dx, qb->ymin,    qb->xmax,    qb->ymin+dy );
}

void qt_node_create_child( qt_node* node /*, float overlap */, int depth )
{
    qt_box ne, nw, se, sw;

    qt_box_split( &(node->box), &ne, &nw, &se, &sw /*, overlap */ );

    node->child[NE] = qt_node_create( &ne );
    node->child[NW] = qt_node_create( &nw );
    node->child[SW] = qt_node_create( &sw );
    node->child[SE] = qt_node_create( &se );
}

int qt_node_has_child( const qt_node* node )
{
    return  (NULL != node->child[NE]);
}

int qt_node_has_data( const qt_node *node )
{
    return  (NULL != node->pconflicts );
}

void qt_search_nodes( qt_node *curr_node, qt_box *sbox, basic_node *results[], int *idx, int max_idx )
{
    if( *idx >= max_idx )  return;

    if( qt_box_overlaps( &(curr_node->box), sbox ) )
    {
        if( qt_node_has_data( curr_node ) )
        {
            results[*idx] = curr_node->pconflicts;
            ++(*idx);
        }

        if ( qt_node_has_child( curr_node ) )
        {
            qt_search_nodes( curr_node->child[NE], sbox, results, idx, max_idx );
            qt_search_nodes( curr_node->child[NW], sbox, results, idx, max_idx );
            qt_search_nodes( curr_node->child[SW], sbox, results, idx, max_idx );
            qt_search_nodes( curr_node->child[SE], sbox, results, idx, max_idx );
        }
    }
}

qt_node* qt_insert_node( quadtree *qt, qt_node *parent, basic_node *node, qt_box *box, int *depth )
{
    if( qt_box_is_inside( box, &(parent->box) ) )
    {
        if( ++(*depth) < qt->depth )
        {
            if( !qt_node_has_child( parent ) )
                qt_node_create_child( parent /*, qt->overlap */, (*depth) );

            if( qt_box_is_inside( box, &(parent->child[NE]->box) ) )
                return  qt_insert_node( qt, parent->child[NE], node, box, depth );

            if( qt_box_is_inside( box, &(parent->child[NW]->box) ) )
                return  qt_insert_node( qt, parent->child[NW], node, box, depth );

            if( qt_box_is_inside( box, &(parent->child[SW]->box) ) )
                return  qt_insert_node( qt, parent->child[SW], node, box, depth );

            if( qt_box_is_inside( box, &(parent->child[SE]->box) ) )
                return  qt_insert_node( qt, parent->child[SE], node, box, depth );
        }

        // inserts into parent since it can't be included in any child nodes */
        node->next = parent->pconflicts;
        // nodes->prev = parent->pconflicts->prev;
        // parent->pconflicts->prev = nodes;
        parent->pconflicts = node;

        return  parent;
    }

    return  NULL;
}

qt_node *qt_insert( quadtree *qt, basic_node *node, qt_box *box )
{
    int depth = -1;

    return  qt_insert_node( qt, qt->root, node, box, &depth );
}

void qt_search( quadtree *qt, qt_box *sbox, basic_node *outlst[], int *outinidx, int inmaxidx )
{
    qt_search_nodes( qt->root, sbox, outlst, outinidx, inmaxidx );
}

/* EOF */
